# R script for (re-)producing maps Figure 1 & 2a
# 2025-05-15

# Make sure the working directory is set to where this script is at
# Or the correct local alternative for the user

rm(list = ls())

# load packages ----

library(sf)
library(ggplot2)
library(dplyr)
library(haven)

# load main data ----

gemeiden_priests <- haven::read_dta('0_georeferenced_hist_gemende.dta') 
priests_repressed <- gemeiden_priests |>
  select(basic_id = id_hist_gemeinde, name_priest)  |>
  filter(!is.na(basic_id)) |>
  filter(name_priest!=".", !is.na(name_priest)) |>
  mutate(name_priest2 = trimws(tolower(name_priest))) |>
  group_by(basic_id) |>
  summarise(
    n_repressed = length(unique(name_priest2))
  ) |>
  ungroup()
gemeinden_all <- gemeiden_priests |>
  select(basic_id = id_hist_gemeinde, name_diocese, lat_4326, long_4326) |>
  distinct() |>
  filter(!is.na(lat_4326), !is.na(long_4326))

bavaria_dta <- left_join(gemeinden_all, priests_repressed, by = 'basic_id') |>
  mutate(is_repressed = ifelse(is.na(n_repressed), 0, 1)) 
bavaria_sf <-  sf::st_as_sf(bavaria_dta, coords = c("long_4326", "lat_4326"), crs = 4326)

gemeinden_sf <- sf::st_read("Shapefiles/Gemeinden_catholic_share_1933/Gemeinden_catholic_share_1933.shp")
gemeinden_sf <- sf::st_transform(gemeinden_sf, crs = st_crs(bavaria_sf))
gemeinden_sf$catholic_1[gemeinden_sf$catholic_1==-9999] <- NA_real_

bayern_sf <- sf::st_read("Shapefiles/Bayern_ex/bayern_ex.shp")
bayern_sf <- sf::st_transform(bayern_sf, crs = st_crs(bavaria_sf))

diocese_sf <- sf::st_read("Shapefiles/Dioceses_poly/Dioceses_poly.shp")
diocese_sf <- sf::st_transform(diocese_sf, crs = st_crs(bavaria_sf))

gedenk_sf <- sf::st_read("Shapefiles/Gemeinden_Gedenk_Residen/Gemeinden_Gedenk_Residen.shp")
gedenk_sf <- sf::st_transform(gedenk_sf, crs = st_crs(bavaria_sf))

target_cities <- c("München", "Nürnberg", "Augsburg", "Regensburg", "Ingolstadt", "Würzburg")
gedenk_cities_sf <- gedenk_sf %>%
  filter(Type == "Stadtkreise", Gemeinden_ %in% target_cities)

resident_sf <- sf::st_read("Shapefiles/residentenliste_bayern/residentenliste_bayern.shp")
resident_sf <- sf::st_transform(resident_sf, crs = st_crs(bavaria_sf))

# plot Figure 1 ----

pdf("Figures/Figure_1.pdf", width = 6, height = 5.5)

ggplot() +
  
  # Background layer: full Bayern map, fills holes
  geom_sf(data = bayern_sf, fill = "darkgrey", color = NA) +
  
  # Gemeinden polygons with fill and borders
  geom_sf(data = gemeinden_sf, aes(fill = catholic_1 * 100), color = "grey40", size = 0.075) +
  
  # Diocese borders
  geom_sf(data = diocese_sf, fill = NA, color = "black", size = 0.7) +
  
  # Priest repression points
  geom_sf(data = bavaria_sf |> dplyr::filter(is_repressed==1),
          aes(shape = "Repressed\nPriests"), color = "black", fill = "#FAE1C3", size = 0.4, stroke = 0.3) +
  
  # Diocese names
  geom_sf_text(data = diocese_sf, aes(label = Name), size = 2.5, color = 'black', fontface = "bold.italic") +
  
  # Color scale for Catholic share
  scale_fill_gradient(
    low = "white",
    high = "red",
    limits = c(0, 100),
    name = "Catholic Share, 1933",
    labels = scales::percent_format(scale = 1, accuracy = 1),
    na.value = "darkgrey",
    guide = guide_colorbar(
      order = 2,
      title.theme = element_text(size = 7.5),
      label.theme = element_text(size = 6)
    )
  ) +
  
  # Manual shape scale to show point legend
  scale_shape_manual(
    name = NULL,
    values = c("Repressed\nPriests" = 21),
    guide = guide_legend(
      override.aes = list(size = 2, stroke = 0.3, color = "black", fill = "#FAE1C3"),
      order = 1, 
      title.theme = element_text(size = 8),
      label.theme = element_text(size = 7.5)
    )
  ) +
  
  # Clean theme
  theme_minimal() +
  labs(x = NULL, y = NULL) +
  theme(
    legend.position = c(0.05, 0.6),
    legend.justification = c("left", "top"),
    legend.box = "vertical",
    legend.spacing.y = unit(0.1, "cm"),
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(0, 0, 0, 0),
    # legend.title = element_text(size = 8),  
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank()
  )

dev.off()


# plot Figure 2 -----------------------------------------------------------

pdf("Figures/Figure_A2.pdf", width = 6, height = 5.5)

ggplot() +
  # 1st layer: Bavaria background
  geom_sf(data = bayern_sf, fill = "grey80", color = NA) +
  
  # 2nd layer: Gemeinden borders (transparent fill)
  geom_sf(data = gedenk_sf, fill = NA, color = "grey40", size = 0.3) +
  
  # 3rd layer: Residents 
  geom_sf(data = resident_sf, alpha = 0.9, size = 0.9) +
  
  # 4th layer: City names
  geom_sf_text(data = gedenk_cities_sf, aes(label = Gemeinden_), color = "white", fontface = "bold.italic", size = 2.75) +
  
  # Theme cleanup
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    legend.position = c(0.05, 0.6),
    legend.justification = c("left", "top"),
    legend.title = element_text(size = 8.5),
    legend.text = element_text(size = 8),
    plot.title = element_blank(),
    plot.subtitle = element_blank(),
    plot.caption = element_blank()
  )

dev.off()







